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Abstract 

Given the high degrees of adaptation to specific microhabitats and restricted-range end- 
emism, Goniurosaurus (Tiger geckos) serves as a unique model to study the complex evo- 
lution in lizards. Using phylogenetic analyses, we estimated the first divergence date of 
Goniurosaurus to the Eocene (~45.3 mya). The diversification within four monophyletic 
species groups began in the mid-Miocene between ~ 13.4 and 7.7 mya and continued to at 
least the early Pleistocene (~2 mya). Their ancestor was predicted to originate somewhere 
in contiguous continental Eastern Asia, whereas the current regions in which each mono- 
phyletic Goniurosaurus species group radiated are respectively their own ancestral regions. 
Together with factors of altitudinal gradient and climate conditions, we reconstructed rel- 
evant niche models of Goniurosaurus including ancestral reconstructions. Consequently, 
low elevations were predicted to be the most probable ancestral state for Goniurosaurus 
and all its groups as well. Both climatic niche conservatism and divergence have shaped 
the extraordinary species richness of allopatric Chinese and Vietnamese tiger geckos. In 
terms of endangerment, Goniurosaurus has been considered one of the most susceptible 
lizard groups under severe human impacts, especially climate change. The assessments 
of their niche evolution can provide a science-based pre-signal of vulnerability, thereby 
improving the efficacy of conservation measures to safeguard species of Goniurosaurus in 
the future. Accordingly, almost all closely related species of Goniurosaurus in China and 
Vietnam were identified with a high rate of niche conservatism, which should be included 
in conservation priorities under potential impacts of climate change. 
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Introduction 


It has been acknowledged that ecological conditions are the major factors that determine 
the distribution of species as well as being the main driver of speciation (Darwin and Wal- 
lace 1858; MacArthur 1984; Orr and Smith 1998; Schluter 2009; Glor and Warren 2011). 
In particular, the largest geographic range in which a species can sustain a stable population 
size, is frequently limited by its physiological tolerances to ecological factors defined by its 
fundamental niche (Grinnell 1917; Kearney and Porter 2004). Together with biotic interac- 
tions, dispersal barriers and chance events, species actually occupy only fractions of their 
fundamental niche, termed realized niche that can explain the high diversity of complex 
range-restricted species, especially in reptile groups distributed in the tropics (Hutchinson 
1957; Sterling et al. 2006; Heaney 2007; Mittelbach et al. 2007; Sober6n 2007; Bohm et al. 
2013; Stein et al. 2014; Steinbauer et al. 2016). In endangerment concerns, such reptiles 
are at high risk of extinction, following Criterion B and D2 in the IUCN Red List relative 
to restricted geographic range (Bohm et al. 2013; IUCN 2021). Against the background of 
increasing anthropogenic impacts, their survival fate is severely imperiled (Thuiller et al. 
2005; Ohlemiiller et al. 2008; Bohm et al. 2013; Cahill et al. 2013; Trew and Maclean 
2021; Cox et al. 2022). Consequently, ecological alternations (e.g., climate change) have 
been globally acknowledged to be a looming factor of extinction, rather than being a driver 
of speciation in the past (Thuiller et al. 2005; Schluter 2009; Lavergne et al. 2013; Bellard 
et al. 2014; Smith et al. 2019). 

Given the potential for niche evolution, the level of adaptability in a species to differ- 
ent environmental conditions can be estimated to assess its chances of survival. Accord- 
ingly, natural populations could be subjected to ecology-based divergent selection across 
the landscape under different abiotic factors (e.g. climate) as non-physical barriers favoring 
reproductive isolation and consequent speciation in new habitats by accumulating genetic 
difference (niche divergence). Consequently, the evolved lineages were detected with dis- 
tinct evolutionary responses from their ancestor and a high level of physiological toler- 
ances to niche alternations (Ahmadzadeh et al. 2013; Krehnwinkel et al. 2015; Rato et al. 
2015; Ahmadi et al. 2021). On the other hand, closely related taxa tend to occupy more 
similar niches than distant phylogenetic congeners and retain ancestral niches over time 
(niche conservatism), when their speciation was mainly related to the allopatric distribu- 
tion patterns by variance. Such species may be more susceptible to climate change due to 
constraints in genetic variation caused by previous niche conservatism limiting the likeli- 
hood of adaptation to changing conditions. Therefore, the looming threat containing novel 
conditions may exceed their narrower physiological tolerance, consequently leading to 
extirpation (Wiens and Graham 2005; Hadly et al. 2009; Wiens et al. 2010; Peterson 2011; 
Lavergne et al. 2013; Pyron et al. 2014; Ahmadi et al. 2021). Understanding the potential 
for niche evolution is particularly important in evaluating the status of endangered species 
and thereby proposing science-based measures of conservation (Richardson and Whittaker 
2010; Wiens et al. 2010; Moritz and Potter 2013). Recently, species distribution models 
(SDM) have been employed to identify patterns and mechanisms of niche evolution, ena- 
bling researchers to test causal relationships between ecological traits and diversification 
processes (Evans et al. 2009; Hof et al. 2010; Salamin et al. 2010; Broennimann et al. 
2012; Lavergne et al. 2013; Rato et al. 2015; Ahmadi et al. 2021). 

Due to the high degrees of adaptation to specific microhabitats and restricted-range 
endemism, tiger geckos belonging to the genus Goniurosaurus Barbour, 1908 can serve 
as a unique model to study the complexity of evolution in lizards (Honda et al. 2014; 
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Liang et al. 2018; Ngo et al. 2021b; Grismer et al. 2021). Recently, systematic issues of 
Goniurosaurus were resolved and ongoing discoveries of new species have been updated 
based on integrative taxonomic approaches using morphology and phylogenetics (Liang 
et al. 2018; Grismer et al. 2021; Ngo et al. 2021b; Zhu et al. 2022). In particular, 25 
tiger geckos have been assigned to one of four monophyletic species groups: the kuroi- 
wae group with six exclusively insular species from the Ryukyu Archipelago, Japan; 
the lichtenfelderi group with five species and the /uii group with nine species (whereof 
G. kadoorieorum is recently considered a synonym of G. luii) distributed throughout 
insular and mainland sites in China and Vietnam; and the yingdeensis group with five 
continental species in southern China (Fig. 1) (Nguyen et al. 2009; Nguyen 2011; Liang 
et al. 2018; Qi et al. 2020; Zhu et al. 2020a, 2020b, 2021; Grismer et al. 2021; Ngo et al. 
2021b, 2022a, b, c). Although most of their distributions are adjacent to each other, 
none of them have been recorded to be in sympatry (Orlov et al. 2008; Yang and Chan 
2015; Zhu et al. 2020a, b, 2021; Ngo et al. 2021b). Oceanic archipelagos and disjunct 
karst ecosystems could serve as physical barriers and/or provide unique microclimate 
zones that might impede dispersal and hence gene flow among the Goniurosaurus popu- 
lations (Clements et al. 2006; Sterling et al. 2006; Heaney 2007). Research presented by 
Liang et al. (2018) indeed underlined the central role of Hainan Island, which formu- 
lated species diversification within the lichtenfelderi group. Using an analysis of ances- 
tral habitat reconstruction, Grismer et al. (2021) shed light on the evolution of habitat 
preference in Goniurosaurus. Accordingly, karstic ecosystems were highly supported as 
the most probable ancestral condition for Goniurosaurus as well as for three of the four 
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Fig. 1 Geographic ranges of the four Goniurosaurus species groups and representative locations of 16 sister 
taxa in China and Vietnam 
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species groups (e.g. kuroiwae, luii and yingdeensis), whereas the lichtenfelderi group 
most likely evolved in non-karstic habitats. 

Goniurosaurus is considered one of the most threatened genera under ongoing severe 
human impacts (e.g. unsustainable exploitation and habitat loss) (Ngo et al. 2019). Fur- 
thermore, as highly endemic ectotherms with limited dispersal capacity, potential fac- 
tors—including small population size, stochastic events and climate change can also push 
Goniurosaurus species to the brink of extinction (Ngo et al. 2019, 2021a, b, 2022a, b). 
To safeguard their populations, legal regulations listing all Tiger gecko species in CITES 
(Convention on International Trade in Endangered Species of Wild Fauna and Flora) 
appendices have recently come into force (Ngo et al. 2019; CITES Notification No. 
2020/068), whereof 17 species have been considered threatened in the IUCN Red List 
(IUCN 2021). In spite of a wealth of studies on systematics, as well as assessments of 
demography and ecology (Liang et al. 2018; Ngo et al. 2019, 2021a, b, 2022a, b; Grismer 
et al. 2021), understanding of the potential of niche evolution in Goniurosaurus remains to 
be investigated in order to identify priorities in species conservation planning. 

Following evolutionary hypotheses, cladogenesis in Goniurosaurus species from East- 
ern Asia was assumed to be strongly related to past orogenic processes that facilitated 
the geographic separation (e.g., disjunct limestone habitats, islands, rivers) which con- 
sequently led to reproductive isolation (Honda et al. 2014; Liang et al. 2018; Ngo et al. 
2021b, 2022b). However, the simple vicariance might play only an indirect role and cannot 
fully explain the mechanism of speciation within Goniurosaurus. Combining a dated phy- 
logenetic tree of all Goniurosaurus species together with their distribution data, altitudi- 
nal gradient and sets of climate conditions, we herein test evolutionary models that may 
explain their extraordinary richness in allopatry and high levels of endemism. In terms of 
species conservation planning, we also intend to identify closely related species of Goni- 
urosaurus preserving climatic conditions over time (niche conservatism), that allows us to 
include them in priorities of conservation protection due to the particular vulnerability to 
climate change in the future. 


Materials and methods 
Genetic data and estimation of divergence times 


Genomic DNA was extracted from muscle tissue samples, using a DNA extraction kit from 
Tiangen Biotech (Beijing) Co., Ltd. Primers used for 16S were r16S-SL (5'’-GGTMMYGC 
CTGCCCAGTG-3’) and l6sbr-H (5'-CCGGTCTGAACTCAGATCACGT-3’) (Palumbi et al. 
1991), for cytb the primers were L14731 (5'-TGGTCTGAAAAACCATTGTTG-3’) (Honda 
et al. 2014) and H15149m (5’-GCMCCTCAGAAKGATATTTGYCCTCA-3') (Chambers 
and MacAvoy 1999), for CMOS the primers were FU-F (5'-TTTGGTTCKGTCTACAAG 
GCTAC-3’) and FU-R (5’-AGGGAACATCCAAAGTCTCCAAT-3’) (Gamble et al. 2008), 
and for RAG1 the primers were R13 (5'-TCTGAATGGAAATTCAAGCTGTT-3’) and R18 
(5’°-GATGCTGCCTCGGTCGGCCACCTTT-3’) (Groth and Barrowclough 1999). The PCR 
procedure was performed with an initial denaturation at 94 °C for 5 min, 35 cycles of 94 °C 
for 30 s, 55 °C for 30 s and 72 °C for 1 min, followed by a final extension at 72 °C for 10 min 
(Liang et al. 2018). PCR products were purified with spin columns and then sequenced with 
forward primers using BigDye Terminator Cycle Sequencing Kit as per the guidelines on an 
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ABI Prism 3730 automated DNA sequencer by Shanghai Majorbio Bio-pharm Technology 
Co., Ltd. 

We constructed Maximum Likelihood (ML), Bayesian Inference (BI), and Bayesian Evo- 
lutionary Analysis by Sampling Trees (BEAST) phylogenetic trees using a concatenated data 
set composed of 3070 base pairs (bp) of the mitochondrial genes, 16 s (633 bp) and cytb 
(1075 bp), and the nuclear genes, CMOS (472 bp) and RAG1 (890 bp), from 103 specimens 
of 24 species of Goniurosaurus with varying degrees of sequence coverage across the sam- 
ples. Concatenation followed the comparison of separate gene trees to confirm there were 
no major discordances. One species, Eublepharis macularius, served as the outgroup (Gris- 
mer 1988; Jonniaux and Kumazawa 2008) to root all trees. See Grismer et al. (2021) for all 
sequence data and GenBank accession numbers. 

A Maximum likelihood (ML) analysis partitioned by gene was implemented using the IQ- 
TREE webserver (Nguyen et al. 2015; Trifinopoulos et al. 2016) preceded by the selection of 
substitution models using TIM2+F+1+G4 for 16 s and cytb and HKY+F for CMOS and 
RAGI. To avoid over parameterization, protein coding genes were not partitioned by codon. 
One-thousand bootstrap pseudo-replicates via the ultrafast bootstrap (UFB: Hoang et al. 2018) 
approximation algorithm were employed, and nodes having UFB values of 95 and above were 
considered strongly supported (Minh et al. 2013). We considered nodes with values of 90-94 
as well-supported. A Bayesian inference (BI) analysis was carried out in MrBayes 3.2.3. 
(Ronquist et al. 2012) on XSEDE using the CIPRES Science Gateway (Cyberinfrastructure 
for Phylogenetic Research; Miller et al. 2010) employing default priors and models of evolu- 
tion that most closely approximated those selected by the BIC and used in the ML analysis. 
Two independent Markov chain Monte Carlo (MCMC) analyses for each data set were per- 
formed—each with four chains, three hot and one cold. The MCMC simulations ran for 100 
million generations. Trees were sampled every 10,000 generations, and the first 10% of the 
trees from each run from each data set were discarded as burn-in. The parameter files from 
both runs were checked in Tracer v1.6 (Rambaut et al. 2014) to ensure that convergence and 
stationarity of all parameters had effective sample sizes (ESS) well-above above 200. Post- 
burn-in sampled trees from each respective run were combined and 50% majority-rule con- 
sensus trees were constructed. Nodes with Bayesian posterior probabilities (BPP) of 0.95 and 
above were considered highly supported. We considered nodes with values of 0.90-0.94 as 
well-supported. 

An input file was constructed in Bayesian Evolutionary Analysis Utility (BEAUti) v. 2.4.6 
using a relaxed lognormal clock with unlinked site models, linked trees and clock models, and 
a Yule prior and run in BEAST on CIPRES (Cyberinfrastructure for Phylogenetic Research; 
Miller et al. 2010). bModelTest was used to numerically integrate over the uncertainty of sub- 
stitution models of each gene while simultaneously estimating phylogeny using MCMC analy- 
ses. MCMC chains were run for 100,000,000 generations and logged every 10,000 genera- 
tions. The BEAST log file was visualized in Tracer v. 1.6.0 (Rambaut et al. 2014) to ensure 
effective sample sizes (ESS) were well-above 200 for all parameters. A Maximum clade cred- 
ibility tree using mean heights at the nodes was generated using TreeAnnotator v.1.8.0 (Ram- 
baut and Drummond 2014) with a burn-in of 1000 trees (10%). Nodes with BPPs of 0.95 and 
above were considered strongly supported. We considered nodes with values of 0.90-0.94 as 
well-supported. 
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Ancestral range estimation 


The resulting BEAST tree was used to estimate the ancestral range at each node using 
the R package BioGeoBEARS (Matzke 2014). We followed the dating scheme of Jon- 
niaux and Kumazawa (2008) in order to compare it with the BEAST and treePL (Smith 
and O’Meara 2012) analyses of divergence times of Liang et al. (2018). BioGeoBE- 
ARS allows for both probabilistic inferences of ancestral geographic ranges and sta- 
tistical comparisons of range expansion from different models in a likelihood frame- 
work employing the Akaike Information Criterion (AIC) to allow the data to choose 
the best fitting model. Available models in BioGeoBEARS include a likelihood version 
of the parsimony-based Dispersal Vicariance Analysis DIVA (“DIVALIKE”) (Ron- 
quist 1997), the likelihood-based Dispersal-Extinction Cladogenesis (DEC) model of 
the LAGRANGE program (Ree and Smith 2008), and the Bayesian-based BayArea 
(““BAYAREALIKE”) (Landis et al. 2013). Additionally, each model incorporates 
founder-effect speciation (+J) which is purported to be particularly important when 
reconstructing biogeographic scenarios of insular lineages (Matzke 2014). For a geogra- 
phy input file, we classified Goniurosaurus species into distinct areas of the four mono- 
phyletic groups’ current geography following the division of Liang et al. (2018), includ- 
ing (1) Guangdong Province, China (C) with four species of yingdeensis group (Qi et al. 
2020); (2) Hainan Island, China (H) with four species of lichtenfelderi group (except 
for G. lichtenfelderi); (3) Ryukyu Archipelago, Japan (R) with six species of kuroi- 
wae group; and (4) Vietnam and China (V) from Red River fault in northern Vietnam, 
northward to southern Guizhou Province, China, with nine species of Juii group and G. 
lichtenfelderi. Each species occurs in only a single region and as such was allowed to 
assign only a single selected area in the analysis (Liang et al. 2018; Grismer et al. 2021; 
Ngo et al. 2021b). 


Occurrence records 


A total of 223 occurrence records of 16 Goniurosaurus species in China and Vietnam 
were collected from our field surveys during the two last decades by ourselves and oth- 
ers (Fig. 1). In order to improve the quality of predictions based on our models, which 
may be affected by geographic bias as results of duplicates representing pseudo repli- 
cates, we randomly selected only one record within each 1 km square by using spatial 
filtering functions in R v 3.1.2 (RStudio Team 2018; Aiello-Lammens et al. 2015). Con- 
sequently, 160 representative records of three Goniurosaurus species groups (including 
lichtenfelderi—44 records, [uii—79 records and yingdeensis—38 records, see details in 
Supplementary information Table $1) were used for further analyses. 


Ancestral elevational range reconstruction 


Based on significant changes of vegetation as Bain and Hurley (2011) suggested for the 
study area, the elevational ranges were divided into three different levels, namely low 
(below 300 m), medium (310-800 m) and high (above 800 m). The coding of eleva- 
tion level for each Goniurosaurus species, which was determined from the literature 
and field observations of the authors, was extracted from coordinate data of occur- 
rences using GPS equipment (Supplementary Information Fig. $1). The classification 
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was mapped onto the BEAST tree of Goniurosaurus using a stochastic character map- 
ping (SCM) implemented in the R package phytools (Revell 2012) in order to derive 
probability estimates of the ancestral states at each node. A transition rate matrix was 
identified that best fit the data by comparing the corrected Akaike Information Criterion 
(AICc) in the R package ape (Paradis and Schliep 2018). Three transition rate mod- 
els were considered: a 2-parameter model having different rates for every transition 
type (the ARD model); a single-parameter model with equal forward and reverse rates 
between states (the symmetrical rates SYM model); and a single rate parameter model 
that assumes equal rates among all transitions (ER). Lastly, an MCMC approach was 
used to sample the most probable 1000 trait histories from the posterior using the make. 
simmap() command and then summarized them using the swummary() command. 


Climatic niche comparisons 


In terms of environmental data, 19 bioclimatic variables for current conditions with a 
spatial resolution of 30 arc-second were obtained from Worldclim (http://www.world 
clim.org/, version 2.0) (Fick and Hijmans 2017). To mitigate the influence of potential 
multicollinearity of predictors on the models, we computed a pairwise Pearson corre- 
lations to identify and eliminate highly correlated climatic variables with r coefficient 
values larger than |0.7 | (Dormann et al. 2013). Our final data set comprised seven 
climatic variables including Mean Diurnal Range (BIO2), Isothermality (BIO3), Max 
Temperature of Warmest Month (BIO5), Mean Temperature of Wettest Quarter (BIO8), 
Precipitation of Wettest Month (BIO13), of Driest Month (BIO14) and Warmest Quar- 
ter (BIO18), which were selected for subsequent analyses. 

To compare the macro-climatic niche among Goniurosaurus groups and sister taxa 
of each group, all pairwise tests were computed by using the package “ecospat” in R 
software (Cola et al. 2017). The niche equivalency and two-way similarity tests, pre- 
sented in detail by Warren et al. (2008) and Broennimann et al. (2012), were calcu- 
lated between pairs of Goniurosaurus groups. The degree of overlap in climatic-niche 
space was further evaluated using Schoener’s D as recommended by Rédder and Engler 
(2011), which varies from 0 (no overlap) to 1 (complete overlap) (Schoener 1970; War- 
ren et al. 2008). However, the niche overlap indices among species-pairs cannot be 
reliably estimated due to the limited availability of records with less than 10 (R6dder 
et al. 2011), which is the case in most Tiger geckos. Following a jackknifing approach 
proposed by Rédder et al. (2011) to overcome the limitation, all available Goniurosau- 
rus records of each monophyletic group were pooled into a dataset and then repeated 
“leave-one-species-out” to calculate Schoener’s D values. We subsequently used the 
value of 1-D as the relative degree of overlap between the omitted one and all remaining 
species. 

Each macro-climatic niche was described in an orthogonal space combining the first 
two principal components in the PCA-env approach (Broennimann et al. 2012). The 
overlap level of given macro-climatic niche spaces was evaluated by direct observations 
on the same environmental background. The proportion of each principal component 
(PC) explaining climatic variance was calculated by the proportion of its eigenvalue to 
the sum. To define the available macro-climatic space, we constructed a minimum con- 
vex polygon (MCP) covering a total of selected Goniurosaurus records using the pack- 
age “adehabitatHR” in R software (Calenge 2006). 


va Springer 


1556 Biodiversity and Conservation (2023) 32:1549-1571 


Species distribution modelling 


For species for which only one or two occurrence records were available (namely G. 
araneus, G. gezhi, G. liboensis, G. kwangsiensis, G. kadoorieorum, G. gollum and G. 
zhoui), we estimated their potential distributions by applying a circular buffer with a 
radius of 10 km around occurrence points with dissolved boundaries in R. In terms of 
SDM, the potential distributions of remaining tiger geckos were predicted by using the 
ensemble of small models (ESMs) combining six modelling techniques: artificial neural 
network (ANN), classification tree analysis (CTA), generalized linear models (GLM), 
generalized additive models (GAM), generalized boosting regression models (GBM), 
and maximum entropy modelling (MAXENT.Phillips). All models were calibrated 
using the “biomod2” package (Breiner et al. 2015; Thuilleret et al. 2016; Cola et al. 
2017). Based on the seven selected climatic variables, we computed a principal com- 
ponent analysis (PCA) in R to derive an orthogonal niche space. Only principal com- 
ponents (PCs) with eigenvalues > 1 were subsequently considered, which were subse- 
quently employed in the SDM (Rato et al. 2015). Although the potential distribution of 
each species was projected within the entire background, all models were only trained 
with presence-only data and 10,000 random points selected within the MCP area. Sub- 
sets of 70% training data and 30% test data were used to model and evaluate the SDMs 
using the two performance indices AUC and TSS (Fielding and Bell 1997; Allouche 
et al. 2006). Values closer to 1.0 indicate better model performance (Breiner et al. 
2015). 


History of niche occupancy 


Following an approach proposed by Evans et al. (2009), the evolutionary history of cli- 
matic niche occupancy of Goniurosaurus species in China and Vietnam was reconstructed 
throughout predicted niche occupancy (PNO) profiles in the “phyloclim” package. The 
suitable probability derived from ESMs of each Goniurosaurus species was re-scaled 
to the sum of 1, and corresponding PC scores binned into 1000 grid cells together with 
designed ones of occurrence-limited species in order to obtain the PNO profile for each 
PC. From these PNO profiles, we resampled 1000 times using the generalized least squares 
method with an assumption of Brownian motion evolution to estimate the climatic spec- 
trum on the BEAST dated tree of Goniurosaurus species in China and Vietnam (Rato et al. 
2015; Ahmadzadeh et al. 2016; Ahmadi et al. 2021). 


Ancestral climatic state estimation 


We tested for phylogenetic signals in climatic niche axes based on climatic PCs and 
our dated phylogenetic tree, using Blomberg’s K (Blomberg et al. 2003) and Moran’s I 
(Miinkemiiller et al. 2012) in the “phylosig” function of the package “phytools” in R. A 
phylogenetic signal was identified with a significance level of p<0.05, indicating that 
Blomberg’s K and Moran’s I are significantly different from zero. In particular, values of 
K range from zero to infinity, where K<1 indicates weak phylogenetic signal, and K> 1 
strong signal, implying that sister taxa are, respectively, more or less divergent in environ- 
mental traits. In terms of Moran’s I, the value reaches to+ 1 indicating a perfect cluster of 
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similar traits, whereas the value of—1 indicates the highest autocorrelation of dissimilar 
traits. 

We illustrated the ancestral state of climatic PCs to assess the pattern of climatic evolu- 
tion of Tiger geckos by using the “contMap” function of the “phytools” package in R soft- 
ware (Revell 2012, 2013). For climatic input data, we averaged the climatic spectrum of 
each climatic PC from the PNO profiles for each Tiger gecko. The best model to describe 
trait evolution was ranked based on the AICs tests for four commonly proposed models: 
Brownian motion (BM), Ornstein—Uhlenbeck (OU), Early Burst (EB) and Lambda, using 
the “‘fitContinuous” function of “geiger” package (Pennell et al. 2014). We computed the 
maximum likelihood (ML) to predict the climatic state following the equation from Felsen- 
stein (1985). 


Results 
Phylogenetic relationships and divergence dating 


The time-calibrated BEAST analysis recovered a phylogeny with well-supported nodes 
(BPP > 90) throughout the tree (Fig. 2). The divergence dates are similar to those calculated 


Me, 
. 
5 


G. splendens 


G. toyamai 


13.3 G. yamashinae 


G. sengokui 


kuroiwae group 


G. orientalis 
G. kuroiwae 
45.3 
® G. zhelongi 
@ G. varius 


G. gollum 


yingdeensis group 
lichtenfelderi group a 
ff / 
Wy 


4 


G. zhoui 


G. kwanghua 


G. lichtenfelderi 
G. hainanensis 
G. bawanglingensis 


G. liboensis 


G. kwangsiensis 


©) Ryukyu Archipelago, Japan 


@ Guangdong, China G. lui = 
@ Hainan Island, China G. kadoorieorum 8 
oils > 

© Vietnam and China G. huuliensis = 
G. catbaensis — 


G. chengzheng 


a 


G. gezhi 


G. yingdeensis | 


G. araneus 


T T T T 1 


40 30 20 10 0 
Millions of years ago 


Fig.2 The ancestral geographical range of Goniurosaurus species reconstructed under BioGeoBEARS 
using DEC +J model on the dated BEAST tree. Colors indicate the current biogeographic unit of each spe- 
cies. Pie charts with colors at the nodes indicate the estimated probability of each area inherited from the 
ancestor after the immediate cladogenetic process 
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by Liang et al. (2018) and fall well within the range of their calculated highest posterior 
density (HPD). 

Based on our dated phylogeny, the first divergence of Goniurosaurus was estimated to 
happen during the Eocene at approximately ~45.3 mya and the radiation continued across 
East Asia up until the Pliocene. Diversification within four monophyletic species groups 
began in the mid-Miocene between ~ 13.4 and 7.7 mya (million years ago) and continued to 
at least the early Pleistocene ~ 2 mya (Fig. 2). 


Ancestral reconstruction of geographical ranges and elevation 


The BioGeoBEARS model indicated that the DEC +J model recovered the best fit to the 
data and most likely to infer the correct ancestral range at each node being that it had the 
lowest AIC-wt score (Supplementary Information Table S2). The range analyses estimated 
that the respective ancestral regions for each group are the same regions in which each 
group currently inhabits. Liang et al. (2018) considered this to be a recent dispersal event. 
However, the geographical origin of Goniurosaurus is not clearly specified. The analysis 
indicated that all the deep nodes outside those of the species groups, recovered an essen- 
tially equal probability that Goniurosaurus originated anywhere across the combined areas. 
Given the young age and relatively recent geographic isolation of the Ryukyu Archipelago 
and Hainan Island from the late Mioence (Sterling et al. 2006; Honda et al. 2014; Liang 
et al. 2018), as for now, we assume that the Goniurosaurus evolved in continental East 
Asia. With additional fieldwork and the potential discovery of new species in the north of 
Guangdong Province, the region of origin could be further refined. 


History of altitudinal transitions 


All Goniurosaurus species have been recorded in areas ranging from 40 to 1000 m elv 
(Supplementary Information Fig. $1). The AIC scores for the three transition models were 
ARD =47.0085, SYM=43.60873, and ER =40.07807. The SCM using the best ER model 
predicted that the “low” level is the most probable ancestral state at the crown tree and all 
group nodes, except for the “high” level at the common node of three lineages (namely G. 
kwanghua, G. hainanensis, G. lichtenfelderi) (Fig. 3). 


Niche overlap comparison 


Climatic niche spaces of three Goniurosaurus groups in China and Vietnam were illus- 
trated in the PCA-env analysis and its first two dimensions accounted for 65.5% of overall 
variance in the climatic data set (whereof, PC1: 38.6%— explained by Mean Diurnal Range 
(BIO2), Precipitation of Driest Month (BIO14) and Warmest Quarter (BIO18); and PC2: 
26.9%—related by Isothermality (BIO3), Max Temperature of Warmest Month (BIOS), 
Mean Temperature of Wettest Quarter (BIO8), Precipitation of Wettest Month (BIO13)) 
(Fig. 4). The analysis revealed no overlaps among their climatic spaces, especially the 
niche space of the yingdeensis group which is entirely separated from the remaining groups 
(Fig. 4). The niche overlap between group-pairs in terms of Schoener’s D only ranges from 
nearly 0 (yingdeensis and two remaining groups) to very limited overlap with D=0.15 
(lichtenfelderi and luii groups). Furthermore, the climatic spaces occupied by the three 
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Fig. 3 Stochastic character map of evolution of elevational occurrence of Goniurosaurus species on the 
BEAST tree. Colored circles at the nodes denote the probability of the ancestral elevation 


groups were neither significantly similar (p>0.05) nor equivalent (p= 1.00) in all group- 
pair comparisons (Supplementary Information Fig. $2, Table $3). 

The Jackknifing approaches omitting records of each species from the pooled data of 
each group documented low “Schoener’s D”—high “1-D” values in some cases, including 
G. lichtenfelderi (1-D =0.998) omitted from the lichtenfelderi group, G. catbaensis (0.844) 
and G. liboensis (0.418) from the Juii group, and G. varius (0.61) from the yingdeensis 
group (Table 1). Our results of the PCA-env analysis further noted in these cases that the 
niche space of each group shrunk significantly (Supplementary Information Fig. $3). 


History of climatic evolution 


The first four climatic PCs from the PCA analysis for SDM account for 90% of overall 
climatic variance, whereof PC1: 39.3% was mainly associated with Mean Temperature 
of Wettest Quarter (BIO8), PC2: 23.9% with Precipitation of Driest Month (BIO14), 
PC3: 18.1% with Max Temperature of Warmest Month (BIOS) and PC4: 7.8% with 
Isothermality (BIO3). These PCs variables were employed together with geographic 
records for the SDMs. The indices of TSS and AUC indicated that the fitted EP-ESMs of 
nine selected Tiger geckos performed well (Supplementary Information Fig. $4). Using 
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Fig. 4 Climatic niche space occupied by Goniurosaurus groups according to the PCA-env analysis includ- 
ing the lichtenfelderi group (granite-stream adapted group)—/uwii group and yingdeensis group (karst 
adapted group), and contribution of climatic variables presented along the first two axes (Correlation cir- 
cle: PC1 and PC2). The solid (100%) and dashed contour (50%) lines illustrate the available macro-climate 
space 


the results of climatic niche models for the potential distribution and circling suitable 
sites of remaining species (Supplementary Information Fig. S5), we reconstructed on 
the BEAST dated tree the ancestral niche occupancy profiles of Chinese—Vietnamese 
Goniurosaurus species, which presented only the niche evolution pattern of divergence 
in given Goniurosaurus groups and intermixed with the conservatism pattern among 
their lineages (Supplementary Information Fig. S5). 

Phylogenetic signals in climatic variables were only detected in PC1l and PC4 
(p-values < 0.05), but these signals were relatively low according to Blomberg’s K 
values, being <1 and Moran’s I values only reached 0.16 (Supplementary Information 
Table S4). The AICc test suggested the Brownian Motion (BM) model was the most 
likely scenario (Supplementary Information Table S4). Consequently, the ancestral state 
estimation also recovered a mix of both climatic patterns (Fig. 6). In particular, the pat- 
tern of divergence was documented among three Goniurosaurus groups in PC4 and 
PC2, and between the yingdeensis group and two remaining groups in PC1 and PC3. 
Given among their lineages, almost Goniurosaurus species followed a conservative 
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Table 1 Niche overlap of 
Schoener’s Index (D) of 16 
Goniurosaurus species of three 
groups in China and Vietnam 


Entire species of each group but D 1-D 
excluding species 


Lichtenfelderi group 


G. bawanglingensis 0.894 0.106 
G. hainanensis 0.959 0.041 
G. lichtenfelderi 0.002 0.998 
G. zhoui 0.975 0.025 
Luii group 
G. araneus 0.861 0.139 
G. catbaensis 0.156 0.844 
G. gezhi 0.94 0.06 
G. huuliensis 0.809 0.191 
G. kadoorieorum 0.930 0.07 
G. kwangsiensis 0.925 0.075 
G. liboensis 0.582 0.418 
G. luii 0.851 0.149 
Yingdeensis group 
G. gollum 0.93 0.07 
G. varius 0.39 0.61 
G. yingdeensis 0.85 0.15 
G. zhelongi 0.65 0.35 


Each calculation excludes one species of each group and 1-D indicates 
the relative degree of overlap between the omitted species and the 
remaining species of each group 


pattern, except for G. lichtenfelderi of the lichtenfelderi group in four PC variables, and 
G. catbaensis in the first three PCs, G. kwangsiensis in PC2 and PC3 and G. liboensis in 
PC1 of luii group (Figs. 5, 6). 


Discussion 
Evolution 


Given estimations of the ancestral range in this study and Liang et al. (2018), the origin of 
the ancestral range of Goniurosaurus was predicted to be somewhere in continental East- 
ern Asia, including the Ryukyu Archipelago and Hainan Island prior to the Eocene, as 
opposed to a specified location. On the other hand, our analyses revealed that the current 
regions in which each monophyletic Goniurosaurus group radiated, are respectively their 
own ancestral region of origin. Dispersals of populations of insular groups (namely kuroi- 
wae and lichtenfelderi) were therefore considered related to prior tectonic events following 
the late Paleocene, such as the collision of India and Eurasia and the clockwise rotation 
of the Philippine Sea plate (Seno et al. 1993; Ota 1998; Teng and Lin 2004; Sterling et al. 
2006; Liang et al. 2018). Subsequent seaways served as a dispersal barrier that prevented 
gene flow among populations. Therefore, vicariant evolution resulted in the formation of 
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Fig.5 Inferred history of the evolution of climatic tolerances in first four climatic PC variables of the PCA 
analysis for species distribution models (SDM) among 16 Goniurosaurus species of three groups (lichten- 
felderi, luii and yingdeensis) in China and Vietnam. Internal nodes represent the mean of climatic toler- 
ances as inferred from the most recent common ancestor of the extant taxa subtended by that node. The 
80% central density of climatic tolerance for each species is indicated by vertical dashed lines and circles 
that mark the respective mean values 


the two insular groups (~late Mioence). The collision also established the land bridge 
accounting for the presence of colonized Indian populations of the Draconinae subfamily 
from Eurasian ancestors during the early to late Eocene (Grismer et al. 2021). 

We used ecological niche models to elucidate the role of the Grinnellian niche (e.g., cli- 
mate conditions) in the diversification of Goniurosaurus in China and Vietnam. Accordingly, 
the niche divergence in all PCs based on climatic variables was strongly correlated with the 
formation of the yindeensis group by the middle Oligocene (Figs. 4, 5, 6). On the other hand, 
the dominant climate of the tropical monsoon toward southeastern China and northern Viet- 
nam was caused by the uplift of Tibetan Plateau as a result of India’s collision with Eura- 
sia as well (Sterling et al. 2006; Miao et al. 2013). The subsequent envelopment of monsoon 
tropical conditions since the Oligocene, might have also driven the cladogenesis of the two 
remaining groups (e.g., lichtenfelderi and luii). In particular, we documented the very low 
overlap between their niche spaces in the PCA-env analysis (Fig. 4) and the discordance in 
their climatic spectrums (e.g., PC2 and PC4) in the PNO estimation (Figs. 5, 6). Furthermore, 
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Fig.6 Ancestral state estimation based on four climatic PCs derived from the PCA analysis for species 
distribution models (SDM) among 16 given Goniurosaurus species in China and Vietnam, reconstructed 
based on the Brownian motion (BM) model, using the contMap function in the phytools R package 
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Grismer et al. (2021) estimated that karst ecosystems are the most probable ancestral habitat 
for most Goniurosaurus groups, except for the lichtenfelderi group which mostly inhabits non- 
karst environments (Grismer et al. 2021). 

Hainan Island is strongly supported as the ancestral range of the lichtenfelderi group (Liang 
et al. 2018). A mixture of karst and granite formations on Hainan Island that may have pro- 
moted the phylogenetic split between G. bawanglingensis and the most recent common ances- 
tor of G. zhoui and G. kwanghua, following the evolution of habitat preference (Liang et al. 
2018; Grismer et al. 2021). The emergence of other species within the lichtenfelderi group 
might generally be associated with climatic cycles during the Pliocene and early Pleistocene 
(Herzschuh et al. 2016). Accordingly, low-to-high transitions along the elevation gradient 
toward highland refugia on Hainan Island occurred during the interglacial mid-Pliocene, and 
thereby tiger geckos (the common ancestor of G. kwanghua, G. hainanensis and G. lichtenfel- 
deri) avoided much warmer climatic conditions. The latest speciation event between G. hain- 
anensis and G. lichtenfelderi in non-karst habitats during the early Pleistocene might relate to 
altitudinal partitioning as well, but vice versa. In particular, the Hainan populations of Goni- 
urosaurus dispersed overwater or overland to occupy low elevations in the Vietnamese main- 
land during sequential glacial periods of sea level decline (Li et al. 2010; Herzschuh et al. 
2016; Liang et al. 2018). Following a pattern of vicariance, these populations were gradually 
isolated by increasing the sea levels, and eventually evolved into G. lichtenfelderi (~2 mya) at 
both island and mainland sites in Vietnam. 

Besides the ocean serving as barrier, the latest speciation event in the lichtenfelderi group 
is also supported by the non-physical barrier of climatic niche limits. We documented the 
pattern of niche divergence in all selected climatic conditions between G. lichtenfelderi and 
Hainan populations of Goniurosaurus (Figs. 5, 6). The role of Grinnellian niche in cladoge- 
netic processes was further assessed in other groups of Goniurosaurus. In general, the low 
phylogenetic signals within Goniurosaurus were strongly influenced by various climate spec- 
tra in the Juii group. In particular, three Goniurosaurus species of the luii group followed the 
pattern of climatic divergence, including G. catbaensis in most climatic conditions (except for 
Isothermality of PC4), G. kwangsiensis in Precipitation of Driest Month (PC2) and Max Tem- 
perature of Warmest Month (PC3), and G. liboensis in Mean Temperature of Wettest Quarter 
(PC1) (Figs. 5, 6). On the other hand, the remaining species in the /uii group and all species 
in the yingdeensis group adapted to their convergent gradients of all climatic conditions, rep- 
resenting a high degree of niche conservatism (Figs. 5, 6). Accordingly, the radiation of their 
common ancestors in Vietnam and China mainland might be enveloped within homogenous 
niche spaces. The subsequent speciation events among these populations of Goniurosaurus 
might be related to potential barriers created by disjunct karst ecosystems, and/ or large riv- 
ers or canyons offering various micro-niches (Clements et al. 2006; Chen et al. 2014; Qi et al. 
2020; Grismer et al. 2021; Ngo et al. 2021b, 2022b; Zhu et al. 2021). For example, the Zuoji- 
ang and Ky Cung rivers might establish an eastern boundary between G. luii and G. araneus 
in China and G. huuliensis in Vietnam, respectively (Ngo et al. 2022b). Given a case of uncer- 
tain taxonomy in Goniurosaurus, the high level of niche conservatism agreed well with the 
morphological and genetic similarities between G. kadoorieorum and G. luii, confirming the 
junior synonym status as proposed by Grismer et al. (2021) and Ngo et al. (202 1b). 


Conservation 


Due to high degrees of endemism and small population sizes, all species of Goniuro- 
Saurus are particularly vulnerable to human impacts (Ngo et al. 2019, 2022b, c; Grismer 


va Springer 


Biodiversity and Conservation (2023) 32:1549-1571 1565 


et al. 2021). In this study, our results suggest that the closely related species within the 
yingdeensis group, the lichtenfelderi group on Hainan Island and some in the /uii group, 
show patterns of niche conservatism and hence appear to be more susceptible to cli- 
mate change (Wiens and Graham. 2005; Hadly et al. 2009; Wiens et al. 2010; Peterson 
2011; Lavergne et al. 2013; Pyron et al. 2014; Ahmadi et al. 2021). Their probability of 
survival in the future has been even reduced decisively by the loss of genetic diversity 
that likely occurs in most Goniurosaurus species due to population declines (e.g., over- 
exploitation) and habitat destruction (e.g., forest fragmentation, quarrying of limestone 
karst and infrastructure) (Frankham et al. 2002; Ngo et al. 2019; 2022b, c). To avoid 
the extinction pressure of climate change, a natural response of species are range shifts 
towards optimal environments (Sinervo et al. 2010). However, the range-restricted tiger 
geckos are unable to migrate larger distances to climatic refugia, which were forecasted 
to be greatest towards higher latitudes in proximity to temperate regions, due to either 
habitat specialization or lacking forest corridors (Williams et al. 2007; Ngo et al. 202 1a, 
2022a). The dense forest canopy, wide elevation gradient and complex karst formations 
likely construct unique micro-refugia containing favorable micro-climate conditions 
which could buffer species against suboptimal climates (Clements et al. 2006; Sterling 
et al. 2006; Ohlemiiller et al. 2008; Suggitt et al. 2018). However, this last survival 
chance is currently getting narrower due to the intensive loss of micro-refugia as the 
result of deforestation and quarrying in karst formations (Clements et al. 2006; Quei- 
roz et al. 2013; CEPF 2020; Ngo et al. 2022b, c). Therefore, higher conservation pri- 
orities and stringent protection should be immediately implemented to prevent these 
niche-conservative species of Goniurosaurus from entering an extinction vortex. How- 
ever, this modelling-based assessment does not necessarily indicate that the remaining 
Goniurosaurus species can overcome environmental changes to adapt. In particular, 
both G. lichtenfelderi and G. catbaensis showing the pattern of niche divergence while 
G. huuliensis following the pattern of niche conservatism, all three of them are poten- 
tially impacted by climate change (Le et al. 2017; Ngo et al. 2021a, 2022a). Although 
approaches to identify the pattern of niche evolution of Goniurosaurus are necessary for 
conservation priorities, the species’ endangered level should be further assessed based 
on comprehensive biological knowledge of population status, ecological characteristics, 
and human impacts. The combination of biological background data is expected to sig- 
nificantly improve the efficacy of conservation plans and actions to safeguard the future 
existence of the species of Goniurosaurus. 
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